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Abstract 

Background: The Orchidaceae is one of the largest families in the plant kingdom and orchid mycorrhizae (OM) are 
indispensable in the life cycle of all orchids under natural conditions. In spite of this, little is known concerning the 
mechanisms underlying orchid- mycorrhizal fungi interactions. Our previous work demonstrated that the non-mycorrhizal 
fungus Umbelopsis nana ZH3A-3 could improve the symbiotic effects of orchid mycorrhizal fungus Epulorhiza repens 
ML01 by co-cultivation with Cymbidium hybridum plantlets. Thus, we investigated the C. hybridum transcript 
profile associated with different beneficial fungi. 

Results: More than 54,993,972 clean reads were obtained from un-normalized cDNA library prepared from 
fungal- and mock- treated Cymbidium roots at four time points using RNA-seq technology. These reads were 
assembled into 16,798 unique transcripts, with a mean length of 1127 bp. A total of 10,971 (65.31%) sequences 
were annotated based on BLASTX results and over ninety percent of which were assigned to plant origin. The 
digital gene expression profiles in Cymbidium root at 15 days post inoculation revealed that 1674, 845 and 
1743 genes were sigificantly regulated in response to ML01, ZH3A-3 and ML01+ ZH3A-3 treatments, 
respectively. Twenty-six genes in different regulation patterns were validated using quantitative RT-PCR. Our 
analysis showed that general defense responses were co- induced by three treatments, including cell wall 
modification, reactive oxygen species detoxification, secondary biosynthesis and hormone balance. Genes 
involved in phosphate transport and root morphogenesis were also detected to be up-regulated collectively. 
Among the OM specifically induced transcripts, genes related to signaling, protein metabolism and processing, defense, 
transport and auxin response were identifed. Aside from these orchid transcripts, some putative fungal genes were 
also identified in symbiotic roots related to plant cell wall degradation, remodeling the fungal cell wall and nutrient 
transport. 

Conclusion: The orchid root transcriptome will facilitate our understanding of orchid - associated biological mechanism. 
The comparative expression profiling revealed that the transcriptional reprogramming by OM symbiosis generally 
overlapped that of arbuscular mycorrhizas and ectomycorrhizas. The molecular basis of OM formation and function will 
improve our knowledge of plant- mycorrhzial fungi interactions, and their effects on plant and fungal growth, development 
and differentiation. 
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Background 

All orchids produce minute, endosperm- lacking seeds 
and are dependent on fungal colonization for germin- 
ation and growth into an underground heterotrophic, 
achlorophyllous stage called a protocorm under natural 
conditions [1,2]. After this fully mycotrophic stage, most 
species shift to autotrophy at adult stage with the devel- 
opment of photosynthetic organs, but approximately 200 
species worldwide stay fully mycoheterotrophic (MH) [3]. 
Recently, some chlorophyllous species in Orchidaceae 
were verified to obtain carbon at adult stage not only from 
their photosynthetic activity, but also from mycorrhizal 
fungi, which were called mixotrophs (MX) [4,5]. Among 
these, several terrestrial orchids such as C. lancifolium 
and C. goeringii were shown to be MX nutritional mode 
through isotopic analyses [6]. The dependency of green 
orchids on mycorrhizal fungi throughout their life cycle 
prompts us to study these mycorrhizal interactions in more 
detail. 

In nature, orchids largely associate with members of 
the 'rhizoctonia' complex, which contains three now taxo- 
nomically disparate Agaricomycetes (/4Hymenomycetes) 
taxa: Sebacinales, Ceratobasidiaceae and Tulasnellaceae. 
Tulasnellaceae are the most frequently found 'rhizoctonias' 
in green orchids [1,3]. Besides orchid mycorrhizal fungi 
(OMF), recent studies also demonstrated that isolates of 
Umbelopsis, as non-mycorrhizal beneficial endophytes, 
probably co-existed in the roots of wild orchids [7,8] . 

C. hybridum collectively contains a large number of 
cultivars with large flowers. Over 20 species in Cymbidium, 
including C. sinense, C. goeringii, C. faberi, C. ensifolium 
and C. kanran, were used in their breeding programs. In 
our previous work, one OMF isolate ML01 and one non- 
mycorrhizal beneficial isolate ZH3A-3 were isolated from 
roots of wild C. spp. adult plants and identified as E. repens 
[ITS accession number of ML01 in GenBank: KJ499806; 
also referred to Tulasnella-like Rhizoctonia] and U. nana 
[ITS accession number in GenBank: HM214447], respect- 
ively. Here, we present that ML01 forms typical mycorrhiza 
with C. hybridum. Individual (ML01, ZH3A-3) or combina- 
tive (ML01 + ZH3A-3, MZ) inoculation significandy im- 
proved the vegetative growth of C. hybridum plantlets 
compared with mock-inoculated controls. 

Although orchid mycorrhizal (OM) symbiosis is an 
important feature in the life cycle of orchids , it has 
long been considered to be an atypical mycorrhizal as- 
sociation, with the fungus deriving little benefit from 
the orchid host [3]. But Cameron et al. [9,10] reported 
that approximately 0.4-3.0% of the labeled carbon 
provided to the orchids was passed to the fungal 
partner [9,10]. Thus adult orchid mycorrhizas may 
represent a truly mutualistic interaction similar to 
ectomycorrhizae (ECMs) and arbuscular mycorrhizae 
(AMs) associations [3]. 



Compared with the numerous transcriptional analyses 
to identify genes involved in AM [11-15] and ECM 
symbiosis [15-17] , the molecular mechanisms behind 
OM symbiosis are still poorly understood [18-21]. Re- 
cently, different gene expression or proteome changes 
of orchid protocorms in response to mycorrhizal fungi 
were studied and mutualistic plant-fungus relationship 
was suggested [19-21]. As for orchids at the adult stage, 
Valadares [20] used RNA-seq to determine the gene ex- 
pression profile of wild Oeceoclades maculata mycorrhizal 
roots in comparison to non-mycorrhizal controls [22]. 
While this study provides new insights into our under- 
standing of the OM symbiosis for adult orchids, the 
symbiosis seems to be highly complex as wild orchids 
generally associate with multiple symbionts even within 
a single root [23]. Although genomic sequence resources 
currently available for orchids are limited, large expressed 
transcripts from Phalaenopsis aphrodite, P. equestris, P. 
bellina, Oncidium Gower Ramsey, Erycina pusilla and C 
sinense were generated in recent years with the rapid 
advances in DNA sequencing technology known as 
next-generation sequencing (NGS or massively parallel 
sequencing) and the development of bioinformatic pro- 
cesses, especially for de novo assembly without refer- 
ence genome information [24-26]. Whereas no genetic 
information on orchid mycorrhizal fungi was available 
when our analyses began, large numbers of assembled 
transcripts and the genome information of T. calospora 
had been recently generated as part of a DOE JGI Com- 
munity Sequencing program. This work has allowed us 
to assign some putative fungal genes to mycorrhizal 
fungi ML01 by Blastn searches of these data. In this 
paper, we used NGS approach and Illumina Digital 
Gene Expression (DGE) technology to obtain genome- 
wide insight into the transcriptomic responses of C. 
hybridum under specific interactions with mycorrhizal 
or/and other beneficial fungi. To our knowledge, this 
represents the first genome-wide gene expression pro- 
file of orchid plants in response to mycorrhizal or non- 
mycorrhizal fungi. The data enabled us to get new 
insight into the mechanisms underlying plant responses 
to mycorrhizal and non-mycorrhizal beneficial fungi. 
The transcriptome for C. hybridum could also be used 
as an important resource for exploring other biological 
processes in orchids. 

Methods 

Biological materials and inoculation 

The micro-propagated plantlets of C. hybridum 'Golden 
Boy' were provided by the Guangdong Key Lab of Orna- 
mental Plant Germplasm Innovation and Utilization, 
China. The propagation procedure for orchid plantlets 
was described in [8]. The two fungal isolates ML01 and 
ZH3A-3 were preserved in the Guangdong Key Laboratory 
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for Innovative Development and Utilization of Forest Plant 
Germplasm. Inocula were produced by growing each strain 
on 90 mm Petri dishes containing potato-dextrose agar 
(PDA) for 10 days. 

The plantlets with 3-4 roots about 2.5-3 cm long 
were then transferred to DE basal medium, an in vitro 
mycorrhizal induction medium in orchids [27]. The co- 
cultivation DE medium contained 1.0 mM CaCl2, 0.5 mM 
MgS0 4 , 1.0 mM K 2 S0 4 , 0.4 mM KH 2 P0 4 , 100 uM FeS0 4) 
25 uM H 3 B0 4) 33 uM MnCl 2 , 2.8 uM ZnS0 4) 1.0 uM 
NaMo0 4 , 140 uM Na 2 EDTA, supplemented with yeast 
extract 1 g • L~ , soluble starch 9 g • IT 1 and agar 6 g • L~ . 
Each plantlet was inoculated with a 7 mm plug excised 
from an edge of an actively growing colony on PDA of 
strain ML01, ZH3A-3 or ML01+ ZH3A-3.The control 
plantlets were mock-inoculated with a plug excised from a 
PDA plate without fungus. All the treated and control 
plantlets were then incubated under a photoperiod of 
16 hrs light at a temperature of 22°C. 

Symbiotic effects on vegetative growth of C. hybridum 

After 45 days post- inoculation (dpi), ninety plantlets of 
each treatment were totally harvested and randomly di- 
vided into three groups, which represented three bio- 
logical replicates. The symbiotic effects on the vegetative 
growth of the hosts were assessed by total increased 
fresh weight and dry biomass of the plantlets. The per- 
centage of root growth (PRG) and mycorrhizal de- 
pendency (MD) were also compared among different 
treatments. MD was calculated by the following for- 
mula: (shoot dry wt. of mycorrhizal plants- shoot dry 
wt. of non-mycorrhizal plants) x 100/shoot dry wt. of 
mycorrhizal plants as described in [28]. For statistical 
analysis of data, SPSS vl3.0 was employed for one-way 
analysis of variance (ANOVA) and comparisons among 
means were made using Duncan's multiple range test, cal- 
culated at p < 0.05. 

Root staining and histology 

The successful inoculation with strains ML01 and ZH3A-3 
and the absence of fungal contaminants in mock-treated 
roots were confirmed by trypan blue staining [8]. To con- 
firm their respective colonization pattern, semi-thin sec- 
tions of root samples for each treatment were prepared 
and stained with a 0.05% aqueous solution of toluidine blue 
for light microscopy (01ympus-BX51, Olympus Corp., 
Tokyo, Japan). 

Reference cDNA library construction 

At 1, 3, 15 and 30 dpi, orchid roots of MZ- treated or 
mock-inoculated control plantlets were harvested by 
carefully removing the substrate and mycelia, and the 
samples were frozen in liquid nitrogen and stored 
at -80°C. Total RNA from each sample was extracted 



separately with Column Plant Total RNAout kit according 
to the manufacturer's protocol (TIANDZ Co., Beijing, 
China). The RNA was treated with DNase (Takara, China) 
and then quantified using a NanoDrop 2000™ spectropho- 
tometer (Thermo Scientific, Waltham, MA, USA). Equal 
amounts of total RNA from MZ- treated and mock- 
inoculated samples at each time point were pooled to pro- 
duce minimum 20 ug total RNA for constructing a cDNA 
library for establishing the reference transcriptome. 

mRNA-sequencing and cDNA de novo assembly 

The library construction and sequencing were performed 
by the Beijing Genomics Institute (BGI) (Shenzhen, 
China) and the general experimental pipeline was sum- 
marized by previously published papers [29,30]. First, 
the cDNA library was pair-end sequenced by Illumina 
HiSeq2000 (90 bp single read module). After removing 
adaptor sequences, empty reads and low quality sequences 
(reads with unknown sequences 'N'), we obtained the 
clean reads for de novo assembly. Three de novo assembly 
software were used including SOAPdenovo with opti- 
mized A:-mer length of 41, ABySS [31] with K-mer length 
of 49 and Trinity [32] (http://trinityrnaseq.sourceforge.net/) 
with default /c-mer length of 25. Sequence directions of the 
resulting unigenes were validated by blastx alignment 
(e-value < 0.00001) between unigenes and public pro- 
tein databases with a priority order of NR (non-redundant 
protein sequences in NCBI), Swiss-Prot, Kyoto Encyclopedia 
of Genes and Genomes database (KEGG), and COG. 
When a unigene happens to be unaligned to none of 
the above databases, ESTScan will be introduced to 
predict its coding regions as well as to decide its se- 
quence direction. 

After sequencing and assembling, we obtained three 
de novo assembled transcriptome database with three 
different assembly methods. The number of assembled 
transcripts, total bases of transcripts, average length of 
all unigenes, N50, number of long-transcripts (>1 kb), 
and number of reads that could be mapped back to tran- 
scripts (RMBT) were compared to evaluate the assembly 
quality with different methods and select the optimum 
reference transcriptome for further study. The expres- 
sion level of unigene was measured by the number of 
clean reads mapped to each annotated unigene normal- 
ized to RPKM (reads per Kb per million reads) [33] and 
adjusted using an additional normalizing factor [34]. 

Functional annotation 

The protein functional annotation of unigenes was 
given according to the known proteins with the highest 
sequence similarity by Blastx to public protein data- 
bases mentioned above. Domain-based alignments 
were carried out against the COG database to predict 
and classify possible functions of unigenes. We further 
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used the Blast2GO program (www.blast2go.org) to get 
GO annotation of unigenes. 

Validation of the reference assembly 

To validate the assembled reference transcripts, 18 
unigenes were randomly chosen for further RT-PCR 
and Sanger sequencing. Gene specific primers were de- 
signed according to the assembled transcripts with pri- 
mer premier software (version 5.0) (Additional file 1: 
Table SI). The reaction was performed in 25 ul, con- 
taining 0.5 (il cDNA template, 12.5 ul Ex Taq Mix 
(TAKARA, Qingdao, China), 1 ul of each 20 uM for- 
ward and reverse gene-specific primers and 10 ul of 
PCR-grade water. Touch-down PCR reactions were 
performed as follows: 95°C for 5 min before cycling 30 
rounds of 94°C for 30 s, 65°C for 30 s (decreasing 0.5°C in 
every cycle) and 72°C for 1 min 30 s - followed by 
15 cycles of 94°C for 30 s, 54°C for 30 s and 72°C for 
1 min 30 s and finally a 72°C extension step for 10 min. 
The specificity of products was inspected by 1% agarose 
gel electrophoresis before they were purified with Ezgene™ 
Gel/PCR Extraction Kit (DC3511-02, Biomiga). The puri- 
fied PCR products were Sanger-sequenced by the 
Beijing Genomics Institute (BGI) (Shenzhen, China). 

Identification of differentially expressed genes (DEGs) 
using lllumina short reads 

To perform quantitative comparisons of the transcript 
variation with respect to different beneficial fungi, total 
RNA samples generated from roots treated with ML01 
(CyEX21), ZH3A-3 (CyEX22) or ML01 + ZH3A-3 (CyEX23 ) 
and mock-inoculated controls (CyEX20) at 15 dpi were 
used to construct un-normalised cDNA libraries which 
were sequenced by lllumina HiSeq2000 (100 bp single 
read module). The cDNA library construction, sequen- 
cing, transcript profiling and mapping back to the 
newly assembled reference transcriptome with Trinity 
were carried out as a custom service (Beijing Novogene 
Bioinformatics Technology Co., Ltd., Beijing, China). 
The expression level of each assembled transcript in 
different samples was measured as RPKM values. 
SeqMap [35] was used for read mapping and rSeq [36] 
was applied for RPKM based expression measurement. 
Abundance data from different samples were collected 
for each transcript. In our work, the changes in rela- 
tive abundance for unigenes between two samples 
were screened with the threshold: q-value < 0.005 and 
log2Ratio > 1 or < -1. 

Quantitative RT-PCR (qRT-PCR) validation 

qRT-PCR was further carried out to analyze the ex- 
pression of twenty- six DEGs identified by RNA-seq 
technology. Gene specific primers were designed ac- 
cording to the cDNAs with Primer Premier software 



(version 5.0) (Additional file 2: Table S2). For RNA ex- 
traction of each treatment, three biological replicates 
were collected independently and immediately frozen 
in liquid nitrogen. One microgram of total RNA was 
reverse-transcribed in a 20 uL reaction mixture from a 
PrimeScript II 1st Strand cDNA Synthesis Kit (Takara, 
Dalian, China). qRT-PCR was performed in a 25 uL re- 
action mixture containing 2 x SYBR Master Premix Ex 
Taq II 12.5 uL (Takara, Dalian, China), 1 uL of cDNA 
template (1:5 dilution), and 1 uL of each correspond- 
ing primer for the gene of interest. qRT-PCR of three 
biological replicates for each sample was performed 
using a LightCycler 480 II System with its relative 
quantification software (ver. 1.2) based on the delta- 
delta-Ct method (Roche). qRT-PCR was performed for 
5 s at 95°C, 10 s at 56°C, and 20 s at 72°C. The cDNA 
samples were standardized to two reference genes: 
B21 (Unigenell515-all) and B24 (Unigenel346-all). 

Results 

Co-cultivation of C. hybridum plantlets with different 
fungi and their respective symbiotic effects on host 
vegetative growth 

When co-cultivated with C. hybridum plantlets, the 
hyphae of the two isolates ML01 and ZH3A-3 began 
to grow at 2 dpi and first contacted the hosts' roots at 
4 dpi. After 15 dpi, the co-cultivation media of plant- 
lets inoculated with ML01 changed from transparent 
to opaque, whereas the actively growing white fluffy 
hyphae of isolate ZH3A-3 fully covered the plate. The 
hyphal growth of isolate ZH3A-3 was retarded by 
the co-inoculated isolate ML01, which showed space 
at the colony margin (Additional file 3: Figure SI). 
After 45 dpi, ML01, ZH3A-3 and MZ treatments all 
significantly improved the total fresh weight and dry 
biomass of C. hybridum plantlets compared to the 
control (Table 1). 

Detection the colonization feature of different beneficial 
fungi on roots of C. hybridum plantlets 

In an initial observation, the hyphae of two isolates 
began to contact with orchids roots at 4 dpi. After that 
time point, we stained the host's roots with each treat- 
ment every two days to monitor the colonization 
process of different fungi and the root staining results 
of each treatment at 6, 10, 15 and 30 dpi were pre- 
sented in Additional file 4: Figure S2. At 6 dpi, the 
first contacted roots were stained to confirm the 
hyphal penetration by isolates ML01 or ZH3A-3. After 
15 dpi, orchid roots were extensively colonized by 
these fungi. At 30 dpi, the symbiotic roots of orchid 
plantlets showed uneven staining owing to the rapid 
growth of roots. To confirm the root staining results 
and characterize the exact colonization pattern of 
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Table 1 Effect of inoculation with Epulorhiza repens isolate ML01 and/or Umbelopsis nana isolate ZH3A-3 on the growth 
parameters of Cymbidium hybridum plantlets at 45 dpi 



Treatment 


PRG (%) 


Leaf FW (g) 


Root FW (g) 


Increased FW (g) 


Leave DW (g) 


Total DW (g) 


MD (%) 


control 


22 


0.76 ± 0.09 b 


0.35 ±0.1 3 d 


0.34 + 0.1 1 d 


0.07 ± 0.01 b 


0.09 ± 0.02 c 


0 


ML01 


533 


1.19±0.14 a 


0.49 ± 0.20 cd 


0.79 ± 0.1 7 bc 


0.13±0.04 a 


0.17±0.04 ab 


46.5 


ZH3A-3 


62.5 


1 .00 ± 0.12 a 


0.71 ±0.13 ab 


0.84 ± 0.1 5 ab 


0.12±0.03 a 


0.17±0.04 ab 


41.7 


ML01 +ZH3A-3 


79.0 


1.11 ± 0.1 2 a 


0.89 ± 0.21 a 


0.99 ± 0.1 4 a 


0.14±0.05 a 


0.20 ± 0.07 a 


50.0 



All the data were mean of three biological replicates, ± standard deviation. (FW: fresh weight; DW: dry weight; MD (%}: percentage of mycorrhizal dependency; 
PRG (%): the elongated root number/total root number x 100). Values followed by different lower-case letters within a column are significantly different at P < 0.05 
according to DMRT. 



these beneficial fungi, we further analyzed the micro- 
scopic structure of the formed symbionts at 15 dpi. 

Root cross sections of C. hybridum plantlets gener- 
ally included root hairs, velamen (2-3 layers), exoder- 
mis (1 layer), cortex (10-13 layers), endodermis (1 layer), 
pericycle (1 layer), vascular cylinder and pith (Figure 1A). 
At 15 dpi, typical pelotons were formed by isolate 
ML01 in the host cortex cells (Figure ID) whereas 
hyphae of isolate ZH3A-3 only colonized the velamen 
cells at regular intervals (Figure IB). Roots co-inoculated 
with isolates ML01 and ZH3A-3 contained newly 
formed and degenerated pelotons, and also typical ag- 
gregated hyphae of isolate ZH3A-3, demonstrating the 
co-infection of two isolates (Figure 1C). None of these 
fungal structures were observed in control root tissue 
(Figure 1A). 



Root transcriptome profile of C. hybridium generated by 
mRNA-sequencing (RNA-seq) 

Up to date, none of genome sequences for Orchidaceae 
was available in public database. RNA-seq, a cost-effective 
and highly efficient next-generation sequencing technol- 
ogy, provides high throughput and is therefore commonly 
used for de novo transcriptome assemblies in non-model 
species, including orchid species [26,37-39]. To analyze 
the transcript variation in C. hybridium colonized by 
beneficial fungi, pooled total RNA from eight samples 
including MZ-treated and mock-inoculated plantlets at 
each time point were used to generate an informative 
reference transcriptome database. Totally, 4,949,457,480 
base pairs raw data were generated using the Illumina 
HiSeq2000, yielding 54,993,972 clean reads that were 
90 bp in length. The raw reads data are available at the 
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Figure 1 Microscopic structure of the symbionts formed by different C. hybridum- beneficial fungi interaction at 15 dpi. Mock- 
inoculated control (A) and inoculated with isolate ZH3A-3 (B), ML01 +ZH3A-3 (C) and ML01 (D). AH, aggregated hyphae of isolate ZH3A-3; Co, 
cortex cells; dH, degenerated hyphae; Ex, exodermis; Mh, fungal hyphae of isolate ML01; Pe, pelotons; Rh, root hair; Va, vascular cylinder; Ve, 
velamen cells. Bar, 50 um 
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NCBI SRA database with the accession number 
SRA051368. 

Three publicly available assemblers as SOAPdenovo, 
ABySS and Trinity were used to de novo assemble 
short-read RNA-seq data into transcripts. The assembly 
software tools as SOAPdenovo and ABySS were chosen 
because of their previous wide application for the simi- 
lar de novo assembly of RNA-seq data sets [31,38,39]. 
Trinity method (release 20110519), one of the newest 
de novo assembly methods, was reported to reconstruct 
of the majority of full-length transcripts in a sample 
from RNA-seq reads directly, across a broad range of 
transcript levels [32]. Recent publications also showed 
that Trinity had a consistently better performance than 
the other assemblers as SOAPdenovo, ABySS, trans-ABySS, 
Oases with single K-mer value, when measured by transcript 
accuracy, integrity and completeness, and sensitivity to 
assemble transcripts from low to high expression levels 
[30]. The outcomes of three assemblers in our study 
were summarized in Table 2. As shown in Table 2, de 
novo assembly using SOAPdenovo yielded a total of 
133,984 unique transcripts (>150 bp), the largest num- 
ber of transcripts. In contrast to the SOAPdenovo and 
ABySS assembly, Trinity assembled the longest mean 
transcript length (1127 bp), the largest N50 (1614 bp), 
the largest number of long-transcripts (>1 kb, 7,170 
unigenes), and the highest percentage of RMBT (40.72%). 
We therefore elected to use the Trinity assembly to retain 
maximum information for further transcript analysis and 
the Trinity assembly is also available in Additional file 5. 

Annotation 

To validate and annotate the assembled unigenes, the 
16,798 unigenes generated by Trinity were subjected to 
BLASTX searches (E-value < le-5) against public protein 
databases. We found 10,915 (64.98%), 8,486 (50.52%), 
5,720 (34.05%) sequences with matches in NR, Swiss-Prot, 
and KEGG databases, respectively. In total, 10,971 (65.31%) 
sequences were annotated and another 6,007 (34.69%) uni- 
genes showed no homology to known sequences in these 
databases. Among the annotated sequences, approximately 



690 transcripts shared the highest similarity to fungal 
or bacterial sequences in NR or Swiss-Prot database. 
The unigenes homologous to known sequences in pub- 
lic databases were further annotated with GO terms 
using Blast2GO. A total of 3,037 (18.07%) unigenes 
were assigned 6,164 GO term annotations in three cat- 
egories: biological process (27.99%), molecular function 
(36.32%), and cellular component (35.69%) (Figure 2). 
Metabolic process, cellular process and response to 
stimulus rank among the first three categories among 
the biological processes. The molecular functions most 
represented were binding, catalytic activity and trans- 
porter activity. Cell, cell part and organelle are the most 
represented categories in cellular components. 

The global transcript abundance could be deduced 
from a non-normalized cDNA library. The transcripts 
highly abundant in the Cymbidium transcriptome are 
listed in Additional file 6: Table S3. There were 5,721 
unigenes mapped into 119 KEGG pathways. Of the 
5,721 unigenes, 1,261 (22.04%) were related to metabolic 
pathways, 545 (9.53%) to biosynthesis of secondary 
metabolites, 531 (9.28%) to spliceosome, 439 (7.67%) 
to plant-pathogen interaction, 154 (2.69%) to protein 
processing in endoplasmic reticulum, and 109 (1.91%) 
to endocytosis (Additional file 7: Table S4). From the 
C. hybridum transcript database, 4,871of the assem- 
bled sequences were assigned to 10,710 COG annota- 
tions, which could be grouped into 25 categories 
(Figure 3). The largest category was "general function 
prediction only" (16.24%), followed by "transcription" 
(10.92%), "replication, recombination and repair" (9.48%), 
"signal transduction mechanisms" (7.63%), "posttransla- 
tional modification, protein turnover, chaperones" (7.18%) 
and "carbohydrate transport and metabolism" (5.73%). 

Validation of assembled transcripts 

To experimentally confirm that the unigenes obtained 
from sequencing and computational analysis were in- 
deed expressed, 18 unigenes were chosen for further 
RT-PCR and Sanger-sequencing. Blastn results showed 
that identities of the selected reference assemblies with 



Table 2 Summary of de novo sequence assembly 



Software 






Length distribution of unigene (nt) 




N50 
(bp) 


Mean 
(bp) 


All 

unigenes 
(nt) 


Length of 
all 

unigenes 
(nt) 


100-300 


300-500 


500-1000 


1000-1500 


1500-2000 


> = 2000 


Trinity 


1,603 


3,401 


4,624 


2,927 


1,811 


2,432 


1614 


1127 


1 6,798 


18,923,761 




9.55% 


20.24% 


27.53% 


1 7.42% 


1 0.78% 


14.48% 










SOAPdenovo 


87826 


28,196 


14,169 


2,692 


734 


367 


349 


329 


133,984 


44,138,252 




65.55% 


21.04% 


10.58% 


2.01% 


0.55% 


0.27% 










ABySS 


No data 


11,303 


8,406 


1,550 


405 


152 


593 


585 


21,816 


12,784,421 






51.80% 


38.50% 


7.10% 


1 .86% 


0.70% 
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COG Function Classification of CyTr-Unigene.fa Sequence 



3 1000 




ABCDEFQH I JKLMNOPQRST JVWYZ 



Function Class 

Figure 3 COG function classification of C. hybridum unigenes. 



A: RNA processing and modification 

B: Chromatin structure and dynamics 

C: Energy production and corwersion 

D: Cell cycle control, cell division, chromosome partitioning 

E: Amino acid transport and metabolism 

F: Nucleotide transport and metabolism 

G: Carbohydrate transport and metabolism 

H: Coenzyme transport and metabolism 

I : Lipid transport and metabolism 

J: Translation, nbosomal structure and biogenesis 

K: Transcription 

L: Replication, recombination and repair 
M: Cell wall membrane envelope biogenesis 
N: Cell molility 

O: Posttranslalional modification, protein turnover, chaperones 

P: Inorganic ion Iransporl and metabolism 

O: Secondary metabolites biosynthesis, Iransporl and calabolism 

R: General function prediction only 

S: Function unknown 

T: Signal transduction mechanisms 

U Intracellular nafficking, secretion, and vesiculai Iransporl 

V: Defense mechanisms 

W: Extracellular structures 

V: Nuclear structure 

Z: Cytoskeleton 
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Sanger-sequencing results were between 96.02-100% 
(see Additional file 1: Table SI), confirming the reliabil- 
ity of assembled reference transcriptome. These uni- 
genes were also included in the next qRT-PCR analysis. 

Global analysis of differential gene expression in the 
symbiotic interaction 

Having generated a C. hybridum root reference tran- 
scriptome, our next goal was to perform comparative 
gene expression profiles in orchid roots inoculated with 
different beneficial fungi. A total of four independent 
un-normalised cDNA libraries were sequenced using the 
Illumina HiSeq2000 100 bp single read module (see 
Additional file 8: Table S5). For each library, the total 
number of counts for each read was determined and all 
the reads were mapped back to the reference transcrip- 
tome (see Additional file 9: Table S6). In this study, 
10.948 million and 14.859 million reads were uniquely 
mapped, representing between 57.02% -76.41% of total 
reads, thus, providing good coverage for the transcript 
profiles. 

The expression level of each assembled transcript 
sequence in different samples was measured through 
RPKM values, and the significant changes in the tran- 
script abundance were identified applying a q-value < 
0.005 and a two-fold change threshold. Based on these 
criteria, at 15 dpi, 1674, 845 orl743 genes were signifi- 
cantly regulated in C. hybridum root inoculated with 
isolates ML01, ZH3A-3 or MZ, respectively. The gen- 
eral differential expressed profiles in the three inter- 
actions are presented in Figure 4 and Additional file 
10: Tables S7-S9. The top 100 abundant genes in sym- 
biotic roots with positive fold change and top 100 
abundant genes in non-symbiotic roots with negative 
fold change in response to symbiotic fungi were pre- 
sented in Additional file 11: Tables S10 and Sll. OM 
formation specifically induced 57% of the top 100 
abundant genes in symbiotic roots with up-regulation 
pattern. Unigene2357_All (18218.59 RPKM), encod- 
ing metallothionein-like protein type 2, which limited 
oxidative damage through binding to and detoxifying 
excess copper and other metals, was the most abun- 
dant gene in OM symbiotic roots with more than 3.6- 
fold up-regulation. Other OM-induced genes in this 
category were largely involved in inhibition of react- 
ive oxygen species (ROS) (glutathione S-transferase 
and superoxide dismutase [Cu-Zn]), protein metabol- 
ism and processing (ribosomal protein, histone H3.2, 
subtilisin-like protease, serine-type carboxypeptidase, 
preprotein translocase subunit SecY), signaling (calmodulin- 
like protein 23, lectin, adhesive/proline-rich protein), cell 
wall modification (extensin class, beta-mannosidase 4, 
beta-tubulin) and cell defense and recue (class III 
chitinase, ATP synthase). Twenty-three genes in this 



CyE X21vsCy EX20 



/ 231 \ 

\ 857 / 42 / 

\ / 544 \ / 

226 \ i ja / 143 



CyEX23vsCyEX20 CyEX22vsCyEX20 

Figure 4 Venn diagram of differentially expressed genes in the C. 
hybridum- beneficial fungi interaction. Mock-inoculated (CyEX20), 
MLOI-treated (CyEX21), ZH3A-3-treated (CyEX22) or ML01 + ZH3A- 
3 -treated (CyEX23). 

v / 

category were co-induced by three treatments, involv- 
ing in signaling (cytokinin-specific binding protein 1, 
annexin, calcium-binding EF hand family protein and 
GTPase), defense (lipoxygenase, GST, class III perox- 
idase, ubiquitin) and general metabolism (isoflavone 
reductase, alcohol dehydrogenase, enolase, lipase). 
However, 63% of the top 100 abundant genes in non- 
symbiotic roots with negative regulation pattern in 
response to inoculants were shared by three treat- 
ments and another 29% of these genes were specifically 
down-regulated by OM fungi. The most representative 
group abundant in non-symbiotic roots and down- reg- 
ulated by three treatments was attributed to genes en- 
coding main components in photosynthesis system, 
whereas a mannose-specific lectin, with 49905.50 RPKM 
in control roots, was specifically down-regulated by OM 
fungi. 

As shown in the Venn diagram, 544 DEGS were 
shared by ML01-, ZH3A-3- and MZ- treatments and 
857 transcripts were collectively regulated in CyEX21 
and CyEX23 but not in CyEX22 which might specifically 
functioned in OM mycorrhizal formation. Individual 
inoculation with non-mycorrhizal fungus ZH3A-3 spe- 
cifically regulated 143 genes. Another 226 DEGS with 
special expression pattern in CyEX23 probably func- 
tioned in the combinative symbiotic process. All the 
DEGs were grouped according to hierarchical cluster 
analysis (Additional file 12: Figure S3) and twenty dif- 
ferent expression sub-clusters based on k-means were 
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identified (Figure 5). The lists of all the DEGs in 20 
sub-clusters were also presented in Additional file 13. 
The genes with the same regulation pattern in three 
treatments were assigned to subcluster 1, 3, 4, 5, 6, 7, 9, 
10, 11, 13, 15, 17 and 19. The genes in subcluster 8, 12, 
18, 20 were specifically regulated by OM formation. 
Genes drastically up-regulated by OM fungi were included 
in subcluster 8, the functions of these genes were related 
to signaling (serine/threonine protein kinase), reactive 
oxygen species homeostasis and stress-related (GST, 
lectin, pathogenesis-related protein 10-3.2), protein 
metabolism and processing (ribosomal protein LlOe, 
subtilisin, serine carboxypeptidase 1, elongation factor, 
amino-peptidase) and transport (ABCG, oligopeptide 
transporter OPT family). 

GO annotation was also used to compare the func- 
tional distribution of DEGS in different treatments 
(Additional file 14: Figures S4-S6). Among the repre- 
sented biological processes, the most represented items 
were "ion transport" for OM formation and "cellular 
biosynthetic process" for individual ZH3A-3 treatment. 
Among the molecular functions, most represented 
DEGs were categorized as "oxidoreductase activity" for 
OM formation and "cellulose synthase activity" for in- 
dividual ZH3A-3 treatment. Among the cellular com- 
ponent, the most represented DEGS were assigned to 
"extracellular region" for OM formation and "plastid" 
for individual ZH3A-3 treatment. 

To identify the biological pathways regulated during 
orchid- beneficial fungi interaction, we assigned the DEGs 
to the reference KEGG pathways (Additional file 15: 
Tables S12-S14). The top 20 pathways with the most 
representation by unique sequences for each treatment 
were included in Figure S7-S9 (see Additional file 16: 
Figures S7-S9). A large number of unigenes related to 
primary and secondary metabolism were collectively 
regulated by three treatments as 'nitrogen metabolism', 
'flavonoid biosynthesis', 'cyanoamino acid metabolism', 
'pentose phosphate pathway', 'starch and sucrose me- 
tabolism', 'alanine, aspartate and glutamate metabolism', 
'stilbenoid, diarylheptanoid and gingerol biosynthesis', 
'porphyrin and chlorophyll metabolism', 'glyoxylate and 
dicarboxylate metabolism', 'Metabolic pathways' and 
'Biosynthesis of secondary metabolites'. Substantial 
numbers of genes/proteins involved in signaling path- 
way, protein turnover, nutrient transport, hormone ac- 
tion, cell rescue and defense were identified in OM 
symbiosis specifically. 

Validation of DEG profiling with qRT-PCR 

Using additional biological samples of ML01-, ZH3A-3- 
and MZ- colonized roots in comparison to mock- 
inoculated control roots, we performed real-time 
RT-PCR to verily the digital expression profiles of 26 




Figure 5 Twenty regulatory patterns of all the DEGs expressed 
in different C. hybridum- beneficial fungi interactions based on 
K-means method. 

■ y 
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genes corresponding to a range of functional categories 
and regulation patterns in the C. hybridum- beneficial 
fungi interactions. The Pearson correlation coefficient 
was calculated by SPSS to assess the correlation be- 
tween different platforms. Overall, the qRT-PCR values 
were highly correlated with the RNA-seq results, 
which confirmed the general expression pattern for all 
the selected genes in three treatments (Table 3 and 
Additional file 17: Figure S10, R 2 = 0.847, correlation 
was significant at the 0.01 level). 

Base on the qRT-PCR results, we further analyzed 
the functional distribution of all the confirmed genes. 
The function of co-regulated genes in three treatments 
included signaling (CBL-interacting serine/threonine- 
protein kinase 11), membrane proliferation (cellulose 
synthase A catalytic subunit 4), scavenging of reactive 
oxygen species (glutathione peroxidase), general stress- 
responsive (laccase-4-like), nutrient transport (inorganic 
phosphate transporter, major facilitator family protein, 
oligopeptide transporter 1-like), protein processing (xylem 
serine proteinase 1 precursor), transcription factor 
(ethylene-responsive element binding factor 4, NAC 
domain-containing protein 43-like, zinc finger CCCH 
domain-containing protein 66-like) and hormone trans- 
port (auxin efflux facilitator PIN1). DEGs derived exclu- 
sively from OM roots were assigned to signaling (LysM 
domain receptor-like kinase 3-like), transport (equili- 
brative nucleoside transporter, ATPase 7, plasma 
membrane-type-like, peptide transporter PTR3-A-like), 
transcription factor (Unigene 2569_A11), cell defense 
(chitinase 4-like) and hormone metabolism (gibberellin 
20-oxidase).Two genes showed special expression 
pattern in response to ZH3A-3 treatment, including 
Unigene 4239_A11 (putative boron transporter 2) and 
Unigene 13219_A11 (putative trans-membrane copper 
transporter). 

Discussion 

Some different effects of mycorrhizal and non- 
mycorrhizal fungi on the growth of C. hybridum plantlets 

This study demonstrated that the vegetative growth 
of C. hybridum plantlets could be improved by orchid 
mycorrhizal fungus ML01 and/or non-mycorrhizal 
isolate ZH3A-3. It is widely accepted that Tulasnella- 
like Rhizoctonia formed typical mycorrhizae with 
orchids and improved their germination and their seed- 
lings' growth [40,41]. Consistent with these studies, 
ML01 formed typical pelotons in host root cortex and 
significantly improved the total fresh weight (FW) 
and dry biomass (DW) of C. hybridum plantlets com- 
pared to the control. However, we observed that the 
roots of ML01- treated plantlets did not elongate ob- 
viously as the co-cultivation proceeded, which might 
be due to rapid acidification of the co-cultivation 



media (our unpublished data). It was thus inferred 
ML01 increased the total dry biomass of orchid plant- 
lets through improving the absorption of mineral 
nutrient. Non-mycorrhizal fungus ZH3A-3 alone 
showed similar promoting effects on total FW and 
DW of Cymbidium plantlets to that of MZ treatment 
at 45 dpi (Table 1), but we also observed that individ- 
ual ZH3A-3- treated orchids initiated less new roots 
than that of dual- inoculated plantlets, and their 
leaves turned into light green after a prolonged co- 
cultivation. Isolate ZH3A-3 always co-existed with 
Tulasnella-like Rhizoctonia in wild orchid roots and 
the hyphae of isolate ZH3A-3 mainly spread in vela- 
men cells and could not form pelotons in the cortex 
cells of Cymbidium roots. It was also revealed that 
the effects of isolate ZH3A-3 and another OMF iso- 
late (CL01) on hosts' uptake of mineral nutrient were 
complementary in our recent report [8]. Thus, we 
designated ZH3A-3 as mycorrhizal-associated fungi. 

Functional analyses of differentially regulated genes by 
OM symbiosis 

The establishment of OM symbiosis involves the mu- 
tual recognition of host and fungi, hypha penetrating 
the cell walls of root hairs or epidermis and then enter- 
ing the cells of the cortical parenchma, forming charac- 
teristic and complex hyphal coils (pelotons) [2]. As in 
other endomycorrhizal systems, the plant cell mem- 
branes in the infected tissues remain intact at all times, 
excluding the hyphae from the cytoplasmic compart- 
ment, and an interfacial matrix is formed between the 
invading hyphae and the invaginating plant plasma- 
lemma [5]. The fungi further spread either as a result 
of repeated colonization or by hyphae penetrating from 
cell to cell in the root cortex, and the cortex cells may 
be uncolonized or contain active pelotons or clumps 
of degenerating hyphae in different proportions. So, 
Rasmussen et al. [5] proposed that there co-existed 
necrotrophic and biotrophic phase in orchid mycorrhiza, 
and the former referred to peloton degradation and unilat- 
eral nutrient transfer by endocytosis, and the latter, mu- 
tual nutritional exchange [5]. A recent labeling study in 
Spiranthes sinensis clearly indicated that in the symbiotic 
protocorm, hyphal degeneration played a significant role 
in the mass transfer of the elements to plant cells, and 
fungal carbon and nitrogen was also detected to be trans- 
ferred from live hyphae through the interface between the 
symbionts [42]. Owing to the non- synchronous process, 
the symbiotic root system often contains fungal structures 
in different developmental stages (contact, penetrating 
hyphae, intracellular hyphae, active and degenerating 
pelotons), and the functions of DEGs for OM symbiosis 
were widely distributed. 
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Table 3 Comparison of the expression profiles of 26 genes by RNA-seq and quantitative real-time RT-PCR 

Gene ID Annotation M/CK (log 2 fold) Z/CK (log 2 fold) MZ/CK (log 2 fold) 







RNA- 
seq 


q RT- 
PCR 


RNA- 
seq 


q RT- 
PCR 


RNA- 
seq 


q RT- 
PCR 


Unigene 1724_AII 


ethylene receptor homolog2 


1.72 


1 .59 ± 0.24 


1.19 


0.73 ± 0.09 


1.87 


1.51 ±0.05 


Unigene 1869_AI 


LysM domain receptor-like kinase 3-like 


4.26 


4.29 ± 0.06 


0.51 


0.98 ±0.1 9 


3.97 


3.85 ± 0.03 


Unigene 2625_AII 


calmodulin binding protein, putative 


1.19 


0.95 ±0.16 


0.87 


0.38 ±0.1 5 


1.30 


0.64 ±0.10 


Unigene 4130_AII 


sugar transporter ERD6-like 7-like 


1.74 


1.70 ±0.21 


0.84 


1.01 ±0.05 


2.03 


1 .63 ± 0.38 


Unigene 4239_AII 


boron transporter 2 


0.96 


0.58 ± 0.25 


-2.17 


-3.07 ± 
0.06 


0.61 


0.27 ± 0.39 


Unigene 6086_AII 


inorganic phosphate transporter, putative 


3.01 


2.81 ±0.34 


2.49 


2.90 ±0.26 


3.40 


3.06 ±0.10 


Unigene 8895_AII 


Xyloglucan endotransglucosylase/hydrolase 


2.03 


0.58 ± 0.66 


2.55 


1 .74 ± 0.49 


1.89 


-0.07 ± 
0.38 


Unigene 
10786_AII 


equilibrative nucleoside transporter 


3.95 


4.71 ±0.11 


-0.21 


0.42 ± 0.08 


3.94 


4.41 ±0.05 


Unigene 
12197_AI1 


ethylene-responsive element binding factor 4 


2.33 


1 .83 ± 0.53 


3.12 


2.72 ±0.39 


3.37 


2.74 ±0.11 


Unigene 
13078_AII 


glutathione peroxidase 


4.20 


2.62 ±0.1 5 


1.67 


1.04 ±0.27 


3.95 


2.45 ± 0.06 


Unigene 
1321 9_AH 


transmembrane copper transporter 


-6.93 


-1.01 ± 
0.03 


8.05 


8.1 5 ±0.44 


6.51 


4.65 ± 0.06 


Unigene 
13639_AII 


ATPase 7, p asma membrane-type-like 


4.24 


4.39 ± 0.09 


0.52 


0.32 ±0.1 1 


3.94 


3.81 ±0.09 


Unigene 
1441 5_A0 


major facilitator family protein 


-2.13 


-1.55 ± 
0.50 


-0.97 


— 1.13 ± 
0.86 


-2.36 


-1.91 ± 
0.50 


Unigene 
15955_AII 


sorbitol transporter 


2.32 


2.86 ± 0.38 


0.87 


1.54 ±0.21 


2.27 


2.20 ± 0.07 


Unigene 
13573_AII 


cellulose synthase A catalytic subunit 4 [UDP-forming]- 

ike 


4.76 


7.26 ±0.91 


4.94 


7.61 ±0.87 


4.67 


6.50 ± 0.89 


Unigene 2569_AII 


transcription factor, putative 


2.18 


2.68 ± 0.30 


0.00 


0.38 ±0.20 


2.00 


2.41 ±0.32 


Unigene 
14999_AII 


NAC domain-containing protein 43-like 


4.19 


3.73 ± 0.36 


4.60 


4.53 ±0.32 


4.32 


3.09 ±0.1 8 


Unigene 
16297_AII 


CBL-interacting serine/threonine-protein kinase 1 1 


1.80 


1 .34 ± 0.35 


1.61 


1.80 ±0.1 8 


2.17 


1.55 ±0.35 


Unigene 
10727_AII 


zinc finger CCCH domain-containing protein 66-like 


2.22 


2.07 ± 0.24 


3.09 


2.82 ± 0.04 


3.53 


2.87 ±0.01 


Unigene 
15958_AII 


chitinase 4-like 


20.38 


13.75 ±0.1 5 


0.00 


-0.04 ± 
0.34 


20.09 


13.40 ±0.21 


Unigene 4429_AII 


laccase-4-like 


2.35 


2.35 ± 0.29 


2.81 


3.38 ±0.06 


1.96 


1.81 ±0.19 


Unigene 
11039_AII 


auxin efflux facilitator PIN1 


3.42 


3.27 ± 0.39 


3.74 


2.84 ±0.07 


3.77 


3.00 ±0.25 


Unigene 
15779_AII 


Xylem serine proteinase 1 precursor, putative 


1.54 


2.09 ± 0.30 


2.40 


3.21 ±0.13 


1.79 


2.32 ± 0.20 


Unigene 2376_AII 


peptide transporter PTR3-A-like 


2.92 


2.38 ± 0.37 


-1.01 


-1.44 ± 
0.45 


2.79 


1.70 ±0.14 


Unigene 2386_AII 


oligopeptide transporter 1-like 


4.51 


4.90 ± 0.49 


1.66 


2.51 ±0.46 


4.57 


4.86 ± 0.49 


Unigene 
13170_AII 


gibberellin 20-oxidase 


1.66 


2.1 3 ±0.05 


0.03 


0.45 ±0.1 6 


1.14 


1 .00 ± 0.27 


Fold changes of selected C. hybridum trasncripts (±standard deviation) were logarithmic transformed which were 


determined by RNA-seq 


and q RT-PCR in root 



tissues treated with isolates ML01 (M) , ZH3A-3 (Z) and ML01+ ZH3A-3 (MZ) at 15 dpi as compared to mock-inoculated control (CK). 



Genes related to signaling pathway or chitin produced by plant interacting microbes and con- 

LysM-domain containing receptor-like kinases (LysM-RLK) firmed to be involved in signaling during nodule and AM 
were identified as receptors of acylated chitin (Nod factors) development [13,43,44]. In this context, it is noteworthy 
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that two genes (Unigene8696_All, Unigene 1869_A11) 
predicted to encode LysM domain receptor-like kinases 
showed evident up-regulation in OM symbiosis. One of 
them, Unigene 1869_A11, encodes a protein with 56% 
identity to LysM type receptor kinase LYS3 from L. 
japonicus (BAI79284). Its induction was validated by 
qRT-PCR, which measured 16-fold up-regulation in the 
mycorrhizal roots. This was consistent with the notion 
that symbiosis signaling seems to occur throughout the 
process of root colonization by a symbiont and occurs 
in parallel to additional signaling mechanisms that re- 
quire the LysM receptor-like kinases [45]. 

Ca 2+ spiking in the nucleus and perinuclear region of 
root hair cells has been documented as one of the earli- 
est cellular responses after the perception of symbionts 
by host plants [46]. Our digital analysis also revealed 
many transcripts associated with Ca 2+ signals, as annexin 
and calcium-binding EF hand family proteins were co- 
regulated by three treatments. But we also detected some 
transcripts in this group were specifically regulated by OM 
fungi, including five genes encoding calcium-binding proteins 
(Unigene4048_All, Unigene6332_All, Unigenell025_All, 
Unigene3685_All, Unigene3072_AU), two genes encoding 
homologues of calcium-dependent protein kinase (CDPK) 
(Unigene6383_All, Unigene8933_AU), one cyclic nucleotide- 
gated ion channel protein (Unigene5311_All), and one gene 
(Unigene6051_All), encoding a homologue of DMI2 nodula- 
tion receptor kinase from M. truncatula. 

Another interesting finding was that two transcripts 
(Unigenel0121_AU and Unigene456_All, up-regulated 
more than 2-fold) encoding putative GRAS family pro- 
tein with equal 44% identity to LjNSPl (nodulation- 
signaling pathway 1 protein) from L. japonicas, and 
one gene (Unigene7236_All, down-regulated more 
than 2-fold), which encoded GRAS family transcription 
factor with 56% identity to NSP2 (nodulation-signaling 
pathway 2 protein) from M. truncatula, were also spe- 
cifically regulated by OM symbiosis. NSP1 and NSP2 
are required for rhizobium-specific gene expression 
whereas only NSP2 might be required for mycorrhizal 
signaling and NSP1 is replaced with required for 
arbuscular mycorrhization 1 (RAMI) to drive expres- 
sion of mycorrhiza-specific genes [45]. 

Genes involved in cellular organization and protein 
processing 

A substantial number of genes related to cell wall 
modification were found to be co-regulated in three 
treatments, which encoded probable xyloglucan endo- 
transglucosylase, beta-D-xylosidase or cellulose syn- 
thase. Two of these were confirmed by qRT-PCR. 
Unigene 8895_A11, encoding a protein with 75% identity to 
the xyloglucan endotransglucosylase (XET) of Medicago 
trancatula, was co-induced more than 4-fold. Another 



gene (Unigene 13573 AU) was strongly up-regulated more 
than 32-fold upon symbiosis (Table 3). The translated pro- 
tein shares 92% identity to cellulose synthase A (CesA) 
catalytic subunit of Vitis vinifera. CesA proteins are part of 
the cellulose synthase complex in higher plants located in 
the plasma membrane [47]. Guether et al. [11] also con- 
firmed that two transcripts of Lotus, LjCell and LjCesA, 
accumulated specifically in the arbusculated cells and func- 
tioned in cell expansion during arbuscule development 
[11]. The strongly up-regulated transcripts of putative 
cellulose synthase might be involved in plant cell wall 
reinforcement during fungal colonization. 

Similar to AM fungi, orchid mycorrhizal fungi pene- 
trate the plant cell wall but not the plasma membrane. 
As the colonization proceeds, the host plasma mem- 
brane invaginates and the infected cells often contain 
numerous mitochondria, ribosomes and profusely devel- 
oped ER. A large number of sequences that were specif- 
ically induced by OM formation were assigned to protein 
metabolism and turnover, membrane dynamics and cell 
wall synthesis. This group of OM specifically induced 
genes in Cymbidium were thought to be homologues as- 
sociated with AM and ectomycorrhizal symbioses [11]. 

Besides the co-regulated DEGS among three treatments, 
OM symbiosis also specifically induced a high number of 
transcripts encoding cell wall-degrading enzymes, includ- 
ing five glucan endo-l,3-beta-glucosidase, one pectinester- 
ase 3 and one xyloglucan endotransglucosylase/hydrolase 
protein. Concomitantly, three extensin- like proteins 
(Unigenell713_All, Unigenel3380_All, Unigenel6702_All), 
and three putative glycosyltransferases (Unigene 1 29 16 A11, 
up-regulated more than 300-fold; Unigenel2916_All and 
Unigene9981_All, up-regulated more than 4-fold) involved 
in cell wall reinforcement were also found in our dataset. 
Overall, this indicates an importance of wall- strength 
maintenance in OM symbiosis. With regard to membrane 
proliferation, two homologues of A. thaliana syntaxin SYP 
132, Unigenel0103_All and Unigene7281_All, were ex- 
clusively induced 16- and 6- fold by mycorrhizal symbi- 
osis. Syntaxin belong to the superfamily of the SNARE 
proteins that play key roles in membrane fusion events 
during vesicle trafficking. In AM symbiosis, a gene 
(Ljwgs_016013.1) encoding a protein with a 60% iden- 
tity to the A. thaliana syntaxin SYP 132 also showed 
seven-fold up-regulation and was suggested to facilitate 
plant/fungal compatibility in mycorrhizal Lotus roots [11]. 

Protein metabolism and processing seems to be a 
major process in mycorrhizal symbiosis [11,48,49]. Most 
drastically up OM- regulated genes in the present study 
were involved in protein synthesis and processing 
(Additional file 11: Tables S10 and Sll). Concerning 
protein synthesis, the strong induction of different 
ribosomal proteins in response to OM fungi provided 
evidence. A number of transcripts associated with protein 
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turnover were also detected in our analysis, including 
three subtilisin-like proteases (Unigene7254_All and 
Unigenel6547_All, up-regulated more than 1000-fold; 
Unigenel0817_All, up-regulated more than 2000-fold), 
three serine carboxypeptidase (Unigene4706_All, up- 
regulated more than 1500- fold; Unigenel2600_All, 
up-regulated more than 300-fold; Unigene5600_All, 
up-regulated more than 6-fold), and also some peptid- 
ase (Unigenel5372_AU, up-regulated more than 320 fold; 
Unigene2916_All, up-regulated more than 40-fold; 
Unigenell542_All, up-regulated more than 4-fold). 
Some transcripts associated with correct protein folding 
(Unigene4797_All, Unigene7897_All) and degradation of 
mis-folded protein (Unigenelll89_AU, Unigenell236_All, 
Unigenel5698_All, Unigene7179_All, Unigene 6288_A11) 
were also identified as specifically up-regulated in OM 
symbiosis. At the same time, a transcript of putative 
serine-type endopeptidase inhibitor, Unigene 1 63 19_ All, ac- 
cumulated more than 40-fold in the mycorrhizal roots 
compared with controls, which might fine-tune protein 
processing. 

Transport-related genes 

Phosphate transfer to the adult orchids via mycorrhizal 
fungi was well documented, just as in other mycorrhizal 
systems. Recently, labeling experiments also directly 
demonstrated that the fully autotrophic orchid Goodyera 
repens acquires carbon, nitrogen and phosphorous from 
its fungal partner [9,10,50]. Notably, G. repens also 
transfers significant amounts of photosynthate (likely 
greater than 3% of its photosynthetic carbon) back to its 
Cemtobasidium mycobiont - the first direct demonstra- 
tion of a net carbon flow from orchid to fungus [9,10]. 
Tulasnella-like Rhizoctonia were also confirmed to im- 
prove the plant growth of C. spp. by enhancing nutrient 
acquisition [8]. 

In our global profiling analysis, 15 genes associated 
with membrane transport were identified as co-regulated 
in three treatments. Among them, 7 genes involved in 
nutrient transport were co-induced in symbiotic 
roots, including two inorganic phosphate transporters 
(Unigene6086_AU, Unigene6124_AU), a plasma membrane 
ATPase 4 (Unigene6958_All), a BCS1 protein precursor- 
like (AAA ATPase chaperone, Unigene 14840_A11) a sorb- 
itol transporter (Unigenel5955_All), an oligopeptide 
transporter 1-like (Unigene2386_All) and a magnesium 
transporter (Unigenelll88_All), whereas 8 genes were 
co-repressed at this time point which encoded two pu- 
tative copper-transporting ATPase, a copper transporter, a 
probable nitrite transporter, a sulfate transporter 3.1, an 
inorganic phosphate transporter 2-1, chloroplastic, a 
probable sodium-coupled neutral amino acid transporter 
6 and a member of major facilitator family. Of these, the 
expression pattern of five genes was further verified by 



qRT-PCR. The first of these, Unigene6086_All, was co- 
induced more than 5-fold, which encodes an inorganic 
phosphate transporter with 81% similarity to inorganic 
phosphate transporter 1-4 of Arabidopsis thaliana and a 
shared 78% identity to inorganic phosphate transporter 
from Lotus japonicus (LjPT3, LjPT4) and M. truncatula 
(MtPT6). Both deduced amino acid sequences of two 
up-regulated inorganic phosphate transporter genes 
(Unigene6086_All, Unigene6124_All) in our digital pro- 
filing have similarity to reported mycorrhiza-specific or 
-induced phosphate transporters [51] (Figure 6). The 
second gene, Unigene2386_All, which was up-regulated 
between 4- and 16- fold in the symbiotic process, en- 
codes a protein with 63% identity to V. vinifera oligo- 
peptide transporter 1-like. Another two co-induced 
genes (Unigenel5955_All, Unigene4130_AU) were in- 
volved in carbohydrate transmembrane transport and fur- 
ther identified by qRT-PCR (Table 3). Unigenel5955_All 
encodes a putative sorbitol transporter with 71% similarity 
to sorbitol transporter of Zea mays and Unigene4130_All 
encodes a transporter with 60% identity to sugar trans- 
porter ERD6-like 7 of A. thaliana. Finally, a member of a 
major facilitator family (Unigenel4415_All) was con- 
firmed to be repressed more than 4-fold in symbiotic 
roots and its translated sequence shared 68% similarity to 
a major facilitator family protein of A. thaliana. In general, 
the regulation of a number of transporter genes supported 
the notion that extensive nutrient exchange was involved 
in symbiosis of the host and beneficial fungi. 

We also identified a substantial proportion of genes 
involved in transport which were specifically regulated 
in mycorrhizal roots. These included three transcripts 
(Unigenell510_All, up-regulated 100 fold; Unigenel3639_All, 
up-regulated 15 fold; Unigene2212_All, up-regulated 32 fold) 
encoding putative plant plasma membrane proton ATPases. 
The amino acid sequence of the encoded protein by 
Unigenell510_All shares 77% identity with the AM- 
specific Mthalin M. truncatula. The regulation of an- 
other transcript (Unigenel3639_All) was validated by 
qRT-PCR (Table 3) and its translated sequence shared 
70% similarity to H + -ATPase 7 from A. thaliana. In 
OM, like in AM and ECM, nutrients from symbiotic 
hyphae separated at the apoplastic space have to go 
through several membrane barriers before being assimi- 
lated by the partner's cells, so the activated plant 
plasma membrane proton ATPase might be a sign of 
accelerated membrane transport. The other two genes 
(Unigenel3137_All, up-regulated 7 fold; Unigene3152_AU, 
up-regulated 2 fold) were also identified as OM-specific 
which encoded a putative copper-exporting ATPase 
and a plasma membrane calcium-transporting ATPase, 
respectively. 

Another group of mycorrhizal-induced genes identified 
was associated with nitrogen transport. Eight genes of 
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Figure 6 Phylogenetic tree for the amino acid sequences of 
plant phosphate transporters (PiTs) described as mycorrhiza- 
specific or -induced transporters, and two putative PiTs of C. 
hybridum (Unigene6086_AII, Unigene6124_AII) in this study 
were highlighted in red color. The dendrogram was generated 
by Mega 4.0 software using ClustlWfor the alignment and the 
neighbor-joining method for the construction of phylogeny. The 
plant Pi-transporters accession numbers could be obtained in [49]. 



peptide transporter (PTR) were found to be exclusively 
up-regulated between two-fold and more than 300-fold 
compared with controls. The strongest up-regulated PTR 
(Unigene4918_All), encoding a homologue of peptide trans- 
porter PTR1 of M. truncatula with 62% identity, showed 
none expression level in the control and CyEX22 samples. 
The second putative PTR (Unigene8412_All) was 300-fold 
up regulated in mycorrhizal roots at 15 dpi and had a 76% 
similarity to a member of oligopeptide transporter family 
(OPT) from P. trichocarpa. Quantitative RT-PCR analysis 
confirmed that one transcript (Unigene2376_All) was up- 
regulated more than 7-fold in OM symbiosis whereas down 
regulated two-fold in ZH3A-3-treated orchid roots. Besides 
these PTRS, two putative amino acid transporters were also 
detected to be induced by mycorrhizal formation. A 
homology search revealed 70% identity between the 
translated amino acid sequence of Unigene3142_All 



and cationic amino acid transporter 5 (CAT5) of A. 
thaliana. The second gene, Unigene404_All, encodes a 
protein with 81% identity to LHT1 (lysine histidine 
transporter 1) of A. thaliana (thale cress), which was 
regulated by diverse factors. Our data supported the 
notion that OMF played an important role in improving 
N acquisition of host plants [9,27] . 

Two putative monosaccharide transporter genes 
(Unigenel2841_All, Unigene4472_All) were also found 
to be induced in OM mycorrhizal roots. The closest 
homologue of translated sequence by Unigenel2841_All 
in A. thaliana is sugar-proton symporter PLT5 (AtPLT5), 
with 71% identity at amino acid level. Another gene, 
Unigene4472_All, encoded a homologue with 70% simi- 
larity to a tonoplast monosaccharide transporter 2 (TMT2) 
in A thaliana and TMT1/2 probably represents a family of 
proton-coupled antiporters capable of high-capacity load- 
ing of glucose and sucrose into the vacuole. Gaude et al. 
[48] also demonstrated that two putative sucrose trans- 
porters and one hexose transporter were induced in non- 
colonized cortical cells of mycorrhizal roots compared to 
the cortical cells of non-mycorrhizal roots, and suggested 
they might be involved in the uptake of carbohydrate 
from the apoplast into cells in the vicinity of hyphae in 
arbuscular mycorrhizal roots or functioned in the 
proton-coupled uptake of carbohydrate into the cyto- 
plasm from the vacuole [48]. This raises a compelling 
question of what is the exact function of these two pu- 
tative monosaccharide transporters in our study. Fi- 
nally, a gene encoding a putative potassium transporter 
(Unigenel4628_All) was also detected to be 6-fold up- 
regulated in mycorrhizal roots and shared 50% identity 
at the amino acid level with the potassium transporter 
1 of A. thaliana. 

In addition to a number of up-regulated genes probably 
involved in nutrient transport, ATP binding cassette (ABC) 
transporter genes, which are unlikely associated with direct 
nutrient transport, were also identified to be regulated in 
mycorrhizal roots. The strongest up-regulated gene of 
ABC transporter family (Unigene5835_All, up-regulated 
nearly 1,000-fold) encoded a protein with 66% identity to 
PDR3 ABC transporter of Oryza sativa, which was regu- 
lated by many abiotic and biotic factors. Recently, two 
plant half-size ATP binding cassette (ABC) transporter 
proteins (STR/STR2) were identified that are essential for 
AM symbiosis [52] and one full-size ABC transporter gene 
MtABCBl (Medtr8g025810) showed strongly increased 
transcript levels in arbuscule-containing (ARB) cells and 
also in adjacent cells [48]. These ABC transporter proteins 
may be involved in the transport of signaling between plant 
cells during AM symbiosis [48]. Here, these identified ABC 
transporter genes showing different expression patterns in 
orchid mycorrhizal roots may be important for OM 
symbiosis. 
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OM-activated transcripts also included components 
associated with vesicle-mediated transport, which might 
involve in the membrane biogenesis for an intracellular 
accommodation of fungal structures. The translated 
sequence of Unigene7513_All (up-regulated more than 
3-fold) shares 70% identity with VAMP7B (vesicle- associ- 
ated membrane protein 7B) in A. lyrata subsp. lyrata. 
Another gene, Unigenell358_All, encoding a putative 
clathrin assembly protein, was up-regulated more than 
8-fold in mycorrhizal roots. 

Defense -related and phytohormone regulation 

Oxidative burst belongs to early events of plant re- 
sponses to compatible and incompatible fungi. The suc- 
cess of beneficial microbes to colonize plant roots depends 
on their ability to manipulate immune responses of the 
host. This was verified by the findings of transcripts as- 
sociated with scavenging of reactive oxygen species 
(ROS), such as glutathione peroxidase, probable gluta- 
thione S-transferase (GSTF1) and germin-like protein 
(Unigene6491_AU, up-regulated more than 37-fold; 
Unigene5007_All, up-regulated more than 25-fold; 
Unigene4304_All, up-regulated more than 4-fold), 
accumulating in the symbiotic roots. One of these, 
Unigenel3078_All, encodes a protein with 77% identity 
to the glutathione peroxidase of Populus trichocarpa 
and was further validated by qRT-PCR. 

Concomitantly, many transcripts encoding general 
'stress- responsive' proteins (e.g. laccase-4-like; probable 
pectinesterase/pectinesterase inhibitor 34 (Unigenell474_AU); 
bibenzyl synthase (Unigenell60_All, Unigenel836_All, 
Unigene4746_All) were identified to be induced in 
these symbiotic process. The enzymes in the pathway 
towards fungicidal phytoalexin synthesis are produced as 
a general feature in orchids. Reinecke [53] have reported 
that the activity of benzyl synthase, catalyzing the biosyn- 
thesis of bibenzyls and 9, 10-dihydrophenanthrenes of 
Bletilla striata rhizomes, increased substantially upon 
fungal infection [53]. Another interesting finding was 
that two transcripts encoding programmed cell death 
protein 4 (PDCD4, Unigene3283_All and Unige- 
nel4957_All) were down regulated collectively in three 
treatments. PDCD4, a suppressor of gene transcription 
and translation, plays a crucial inhibitory role in several 
types of human tumors through inhibiting autophagy in 
multiple cell types both in vitro and in vivo [54]. Down- 
regulation of these homologues of PDCD4 may pro- 
mote the symbiotic process through suppressing plant 
defense. Furthermore, four transcripts (Unigene4420_AU, 
Unigene8976_AU, Unigene8566_All, Unigenel3020_All), 
encoding homologues to key molecular regulators in 
two-component systems (TCSs) which play important 
roles in responses to environmental stress stimuli [55], 
were found to be down regulated regardless of the 



inoculated fungi species, highlighting the shared strat- 
egy by different beneficial fungi in the symbiosis. 

Apart from these shared DEGs, three chitinase genes 
(Unigenel5958_All, Unigenel6629_All, Unigene8905_All) 
were exclusively induced by OM fungi and the specific ex- 
pression pattern of Unigene 15958_A11 was confirmed by 
qRT-PCR. The translated sequence of Unigene 15958_A11 
shared 59% identity with the chitinase 4-like of Fragaria 
vesca. Chitinases hydrolyze [3-1,4-glycosidic bonds and 
can be induced in pathogenic interactions. But the ex- 
pression of Unigene 15958_AU was not induced by a 
pathogenic Fusarium isolate (our unpublished data). 
Unigene8905_All encoded a protein with 55% identity 
to class III chitinase of M. truncatula and class III 
chitinases were found to be expressed predominantly in 
the arbuscule-containing cells [14,56], where they may 
suppress plant defense by the reduction of chitin-like 
elicitors during the formation of functional symbiotic 
interfaces [57]. The similar role of class III chitinase 
SvChit3 in OM was recently suggested by Perotto [19]. 
Thus, these up-regulated chitinases are probably in- 
volved in the lysis of fungal cell wall and balancing the 
symbiosis. 

Another notable OM-regulated group was mannose- 
specific lectins (Unigene4318_All, Unigene529_All, 
Unigenel4222_All). Among these, Unigene4318_All 
and Unigene529_All showed high-level expression levels 
in mycorrhizal roots with little expression in the un- 
inoculated control and ZH3A-3 treatments. Antifungal ac- 
tivity against a range of phytopathogenic fungi has been 
demonstrated for some of mannose-specific lectins both 
in vitro and in transgenic plants, and they were thus pro- 
posed to be part of defense mechanisms against the invad- 
ing fungus. The strong up-regulation of SvLect3 and 
SvLectS in response to mycorrhizal fungi T. calospora in S. 
vomeracea protocorm was recently confirmed by qRT- 
PCR and this was suggested to be part of a mechanism 
to limit fungal growth in the mycorrhizal protocorm 
[19]. However, another mannose-specific lectin gene, 
Unigenel4222_All, with 49905.50 RPKM in control 
roots, was specifically down-regulated by OM fungi. 
Owing to the multi-functions of plant lectins, the role 
of these lectins in the orchid mycorrhizal process re- 
mains to be established. 

Among plant hormones, ethylene (ET), salicylic acid 
(SA), abscisic acid (ABA), and jasmonic acid (JA) are 
known to be key elements in fine-tuning the plant 
defense response during interaction with other organ- 
isms [58,59]. The inhibitory effect of exogenous ethylene 
in AM development has been documented and the 
colonization of host tissues by mutualistic AM fungi such 
as Rhizophagus irregularis (formerly Glomus intraradices), 
ectomycorrhizal fungus Laccaria bicolor, or endophytic 
fungus Piriformospora indica was inhibited by ET [60]. 
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However, ET signaling is still essential for nodule for- 
mation in Sesbania rostrata and in the initial stages of 
P. indica colonization [61]. Because ethylene-responsive 
transcription factors (ERF) can function as transcrip- 
tional activators and repressors, they are candidates for 
establishing a balanced defense response to the fungus 
without preventing growth and development in AM 
symbiosis [62]. Valadares [20] also suggested the role of 
ethylene in OM regulation when compared the gene ex- 
pression pattern of O. maculata mycorrhizal roots in 
comparison to non-mycorrhizal controls [22]. In our glo- 
bal expression profiling, eight transcription factors (TFs) 
related to ethylene signaling (Table 4) were found to accu- 
mulate in all the symbiotic roots and four genes among 
these were specifically induced in OM roots. Among 
these, the induction of two transcripts (Unigenel724_All, 
Unigenel2197_All) was further validated by by qRT-PCR. 
The translated sequence of Unigenel724_All showed 63% 
similarity to the ethylene receptor homolog 2 (ETR2) of 
Zea mays with two-fold induction. Unigenel2197_All, 
encoding a protein with 45% similarity to ethylene- 
responsive element binding factor 4 (ERF4) of G. max, 
was induced more than 5- fold in the symbiotic roots. As 
with Plett et al. [60], the induction of these TFs is probably 
an adaptive response by the plant so that its growth and 
vigor are not compromised by the fungus in the late stage 
of symbiosis [60]. 

ABA plays a crucial role in biotic and abiotic stress 
responses. Analysis of AM colonization in the sitiens 
ABA-deficient tomato mutant showed that ABA was 
necessary in order to complete arbuscule formation and 
its functionality and to promote sustained colonization 
of the plant root [63]. Subsequent work supported the no- 
tion that ABA deficiency negatively regulates mycorrhizal 
intensity/hyphal colonization directly or indirectiy through 
increasing ethylene levels [64-66]. Accumulation of the 
transcripts related to ABA de novo biosynthesis was 
documented in AM roots [49] and OM green protocorms 
[20]. At the same time, genes associated with ABA catabolism 
also showed increased transcript levels in AM roots [49] and 
OM roots [22]. In our study, 4 transcripts related to ABA 
biosynthesis: 9-cis-epoxycarotenoid dioxygenase (NCED) 
(Unigenel444_All, Unigene2412_All) and zeaxanthin 
epoxidase (ZXE) (Unigenell284_All, Unigenel3062_All), 
were co-suppressed in symbiotic roots. ABA deficiency 
might be involved in restriction of fungal spread in orchid 
symbiotic roots. 

Jasmonic acid has long been implicated in the plant's 
systemic response to pathogen attack [67]. While experi- 
mental results sometimes appear contradictory, most 
authors agree that jasmonate signalling is important for 
AM colonization and development [68]. Plett [60] re- 
vealed that in the interaction between L. bicolor and 
poplar, ethylene and jasmonic acid showed extensive 



transcriptional cross-talk and functioned as negative 
modulators during mutualistic symbiosis [60]. Jasmonates 
belong to a diverse class of lipid metabolites known as 
oxylipins and the key genes in the biosynthesis of oxylipins 
were commonly induced in mycorrhizal roots [59]. Con- 
sistent with these studies, the induction of 6 lipoxygenase 
(LOXs) genes in our study suggest the involvement of the 
oxylipin pathway in regulation of OM symbiosis. In plants, 
there are two main branches of the oxylipin pathway, de- 
termined by two different types of lipoxygenases (LOXs), 
9-LOX and 13-LOXs. The 13-LOX pathway leads to the 
biosynthesis of JA and derivatives, and the relevance of 
the 9-LOX pathway in plant interactions with nematodes 
and pathogens has also been demonstrated recently 
[69,70] . Further research is required to elucidate the exact 
role of these up-regulated LOXs in long-term balance of 
the OM symbiosis. 

Apart from the potential roles for ET, ABA and JA in 
OM symbiosis, three transcripts associated with gibberellin 
(GAs) biosynthesis were also increased in abundance in 
orchid mycorrhizal roots: two genes (Unigenell581_AU, 
Unigene3993_AU) encoding gibberellin 3-oxidase (GA3ox) 
and another (Unigenel3170_All) encoding GA 20-oxidase 
(GA20ox). The translated sequence of Unigenel3170_All 
shares 77% identity to gibberellin 20-oxidase of Triticum 
aestivum and its induction was confirmed by qRT-PCR. 
The product of GA3ox and GA20ox activates GA, suggest- 
ing that GA synthesis is essential for OM formation. In 
barley- P. indica interaction, Schafer [71] also observed 
that two barley mutants impaired in GA synthesis or per- 
ception, reduced colonization by P. indica. Similarly, 
homologues of genes involved in the GA biosynthesis were 
already identified as mycorrhizal-responsive in tomato 
plants at the mature (50 dpi) phase of the interaction [66]. 
Foo [68] reported that biologically active gibberellins sup- 
press arbuscule formation in pea roots, and gibberellin sig- 
naling DELLA proteins are essential for this response. 
Owing to extensive transcriptional cross-talk between GA 
signaling and SA/JA responses, the GA-regulated defense 
response might be critical for establishing these long-term 
mutualisms. 



Transcription factors (TF) 

As detailed above, three treatments initiated common 
stress responses for symbiotic roots, so the majority 
of co-regulated TFs encoded AP2 domain TFs and 
ethylene-responsive TFs (ERF), WRKY TFs, Myb- family 
TFs, basic helix-loop-helix (bHLH) TFs and zinc-finger 
TFs (Table 4). DNA-binding WRKY, AP2-domain and 
ERF TFs are generally involved during pathogen-related 
stress responses [48,61]. The co-regulated TFs may also 
be involved in root developmental process, as are MYB 
family, ERF and NAC family [72]. 
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Table 4 Transcription factors among beneficial fungi- induced genes 


TF family 
or domain 


Induced 


Repressed 


Number 
of genes 


AP2-ERF 


Unigene627_AII*, Unigene6402_AII*, Unigenel 21 97_AII*0, Unigene103_AII*, 
Unigene1545_AII, Unigene10693_AII, Unigene13239_AII, Unigene13630_AII 




4* (4) 


WRKY 


Unigene2569_AII 0 (WRKY19), Unigenel 3 120_AII*(WRKY40), Unigene8345_AII* 
(WRKY5), Unigene8346_Alr(WRKY5), Unigene3004_AII (WRKY39), 
Unigene451 1_AII*(WRKY43), Unigene1640_AII(WRKYl7) 


Unigenel 887_AII (WRKY23) 


5* (3) 


Myb 


Unigenel 0890_AII*, Unigene13524_AII*, Unigene14312_AII* Unigene2947_AII 


Unigenel 2533_AII, Unigenel 61 61 _AII*, 
Unigenel 441 All*, Unigenel 4070_AII* 


6* (2) 


bHLH 


Unigene4212_AII*, Unigene2945_AII* Unigenel 3683_AII*, Unigene14108_AII* 


Unigenel 570_AI1* Unigenel 5208_AII*, 
Unigenel 2890_AII, 


6* (1) 


HOX 


Unigenel 1426_AII*, Unigene888_AII, 


Unigene9904_AII, Unigene7407_AII*, 
Unigene14238_AII 


2* (3) 


NAC 


Unigenel 4999_AII*ffi, Unigene8664_AII*, Unigene16437_AII*, Unigenel 1336_AH* 


Unigene5046_AII* 


5* 


ARF 


Unigene! 0085_AII (ARF9), Unigene9466_AII (ARF9) 




(2) 


GRAS 


Unigene2120_AII(SC)*, Unigenel 01 21_AII, Unigene4453_AII (scarecrow-like), 

I ninQnafiylQI All 1 ninQnolSin All /<Zf~} 1 l-linonol ^^^fl All 

unigeneo4oz_Aii pLj, unigene IjjIz_aii pLj, unigene jcoU_aii, 
Unigene456_AII 




1*(6) 


BTB/POZ 


Unigene2160_AII, Unigene351_AII, Unigene352_AII 


Unigenel 1206_AI 


(4) 


C3H 


Unigene156_AII, Unigene209_AII*, Unigenel 0727_AII*® 


Unigene16285_AII, Unigene81_AII 


2* (3) 


GATA 


Unigenel 952_AII* Unigenel 01 79_AII* 


Unigene1691_AII, Unigenel 021 9_AII* 


3*0) 


C2H2 


Unigene2954_AII*, 


Unigene2352_AII* 


2* 


dof 


Unigene3028_AII*, 




1* 


A20AN1 


Unigene2954_AII*, Unigenel 0730_AII 




1*0) 


CO-like 




Unigene4931_AII 


(1) 


Hsf 


Unigenel 3648_AII*, Unigene4980_AII* 




2* 


Trihelix 




Unigene2892_AII 


0) 


ASG 




Unigene3840_AII 


0) 


WD 


Unigene4225_AII 




0) 


MADs 


Unigene4681_AII 




0) 


MACPF 


Unigene8197_AII 




0) 


BLH 


Unigene8238_AII 




0) 


A-2b 


Unigenel 3648_AII*, Unigene4980_AII 




1*0) 


UNE 


Unigene4910_AII 


Unigenel 0862_AII* 


1*0) 


rf 


Unigene14563_AII* 




1* 


MYC 


Unigenel 5000_AII* 




1* 


RAX 


Unigenel 589_AII* 




1* 


NF-YC 


Unigene16018_AII* 




1* 


bZIP 


Unigene4820_AII, 




(1) 


Lists of TFs i 
parenthesis 


n this table included specifically regulated TFs in OM formation and * co-regulated genes 
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We also detected GRAS family and ARF TFs specif- 
ically induced by OM formation. Four transcripts 
(Unigene4453_All, Unigene6482_All, Unigenel3312_All, 
Unigenel5630_All) encoding for GRAS members were 
putative SCARECROW (SCR) protein encoding genes, 
considered to be genetic determinants of root identity. 
Previous studies also confirmed that members of the 
GRAS family are essential for nodule development [11] 
and they may be important for the regulation of gene 



expression in AM roots of M. truncatula, L. japonicus and 
5. lycopersicum [11,12,49]. Two putative ARF transcripts 
were also identified to be up- regulated in orchid mycor- 
rhizal roots. The induction of putative ARFs has also been 
found during AM symbiosis in maize, rice and M. trunca- 
tula, but not in L. japonicus [72]. And in AM interaction, 
mycorrhizal roots exhibit morphological characteristics 
such as an increase in the number of lateral/fine roots 
during early growth phases, similar to auxin-treated roots 
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[49]. Recently, Sukumar (2012) and Fusconi [73] summa- 
rized a compelling supportive role for involvement of 
auxin pathways in beneficial plant-microbe interactions 
[73,74]. Chutima [75] reported that orchid-associated 
fungi Tulasnella sp. produced high levels of indole-3- 
acetic acid (IAA) with all biological activities in culture 
medium supplemented with 2 mg/ml of L-tryptophan 
[75]. The IAA production by isolates ML01 was also iden- 
tified in our research group (unpublished data). Further- 
more, Unigenell039_All, predicted to translate a protein 
with 66% similarity to an auxin efflux facilitator S1PIN1 of 
Solarium lycopersicum, showed 15-fold up-regulation in 
the symbiotic roots. All these data point to the probable in- 
volvement of auxin signaling pathway in OM functionality. 

DEGs specifically regulated by non-mycorrhizal fungus 
ZH3A-3 

The hyphae of non-mycorrhizal fungus ZH3A-3 colo- 
nized the velamen cells of orchid roots and individual 
inoculation significantly improved the root length and 
FW of C. hybridum plantlets. As detailed above, 845 
genes were significantly regulated in C. hybridum roots 
inoculated with ZH3A-3 alone and most of these DEGs 
were shared with those of ML01 and/or MZ treatment. 
Apart from these shared DEGs, we also identified 143 
genes were specifically regulated by ZH3A-3 treatment. 

The largest number of up-regulated transcripts by 
ZH3A-3 treatment was associated with orchid root de- 
velopment. These genes are mainly related to hormonal 
balance and transcription regulation. Among plant 
phytohormones, auxin, cytokinin, and ethylene play 
important role for root development. Consistent with 
this, we detected the up-regulation of many corre- 
sponding genes, which encoded one cytokinin dehydro- 
genase (Unigenel4465_All), two ERF (Unigenel5166_All, 
Unigene7827_All), and three ABC transporter B fam- 
ily members (Unigenel6320_All, Unigene5921_All, 
Unigenel4693_All). Auxin plays a pivotal role in many 
plant- beneficial microbe interactions [74], and it is well 
established that higher auxin/cytokinin ratio determines 
the root initiation. A recent study also clearly demon- 
strated that OsCKX4, a cytokinin oxidase/dehydrogenase 
(CKX) family gene, integrates cytokinin and auxin signal- 
ing to control rice crown root formation positively [76]. 
So, the up-regulation of Unigenel4465_All, which en- 
codes a homolog with 69% identity to CKX of A. thaliana, 
may suggest its positive role in orchid- ZH3A-3 inter- 
action. Polar auxin transport (PAT) is essential for estab- 
lishment of auxin gradient that spatially defines the 
quiescent center (QC) and consequently the fate of 
neighboring stem [77]. ABC transporter B family mem- 
bers were probably involved in auxin acropetal or ba- 
sipetal transport [78], thus, the up-regulation of three 
ABC transporter B family members appears to converge 



on the auxin pathway to modulate root development. 
As shown in Table 4, the putative MYB-TFs ranked the 
third largest number of the regulated TFs in the 
mycorrhizal symbiosis, and most were shared by three 
treatments, which suggested their essential role in re- 
sponse to symbiotic fungi. Recently, a convincing re- 
search supported that a putative MYB-like TF of L. 
japonicus {LjMAMI) which was AM-responsive may 
also have non-symbiotic functions: i.e. root growth 
[79]. In our study, two transcripts (Unigene2947_All 
and Unigenel589_All), encoding homologues of Myb 
family TFs of O. sativa Japonica Group and A. thaliana 
with 60%- 80% identity, were also detected to be up- 
regulated 4-fold by ZH3A-3.The translated sequence 
of Unigene2947_All also shared 46% similarity to 
LjMAMI at amino acid level. Another up-regulated TF 
related to root development was Unigene489_All (up- 
regulated 2-fold), which encodes a protein with 45% 
similarity to SCARECROW-like protein of Z. mays. 

We also identified a substantial number of genes re- 
lated to general disease resistance processes which were 
specifically up-regulated by ZH3A-3 treatment. These 
included 3 transcripts related to signaling (leucine-rich 
repeat transmembrane protein kinase), 5 transcripts 
associated with balance of reactive oxygen (respiratory 
burst oxidase homolog, class III peroxidase, L-ascorbate 
oxidase, monocopper oxidase-like protein), and 12 tran- 
scripts involved in secondary metabolism (phenylalanine 
ammonia-lyase, DnaJ homolog, disease resistance re- 
sponse protein, isoflavone 2'-hydroxylase, syringomycin 
biosynthesis enzyme). The increased activity of chitinases 
were generally found in nonspecific, broad-spectrum de- 
fenses, but we did not detect the transcript accumulation 
of chitinase genes in CyEX22.The fungal cell wall of iso- 
late ZH3A-3 always ruptured as the cultivation proceeded 
(our unpublished data), which might be the reason. The 
activation of basal defense mechanisms of orchid plants 
by beneficial fungus ZH3A-3 did not affect the symbiotic 
effects, as described in [80] . 

As ZH3A-3 colonization induced dramatic structural 
changes in root cells and promoted the root growth of 
orchids, transcripts accumulated in the symbiotic roots 
were also related to carbohydrate and amino acid metab- 
olism, protein synthesis and processing, and structural re- 
organization. Furthermore, a transcript (Unigenel5068_All) 
involved in water transport was also up-regulated, which 
encodes a homolog of aquaporin PIP1-2 of Z. mays with 
87% identity. 

Taken together, we confirmed the findings which 
suggested that the development and metabolism of plant 
symbionts are generally driven by differential regulation 
of transcriptional regulators, signal transduction, and 
metabolic pathways, rather than by expression of 
symbiosis-specific genes [81]. 
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Detection of putative fungal genes in symbiotic roots 

Owing to the distinct phylogenetic relationship between 
isolate ML01 and ZH3A-3, we could detect 163 possible 
fungal genes combining the annotation information of 
de novo assembled transcripts with the comparative 
expressed level in different symbiotic roots (Additional 
file 18: Tables S15 and S16). Of these 163 genes, 53 were 
probably genes from isolate ZH3A-3, and the rest, from 
ML01. The origin of putative ML01 genes was further 
confirmed by Blastn search against the T. calospora 
draft genome or T. calospora mycelium ESTs with an 
E value > le~ 20 (http://genome.jgi.doe.gov/pages/blast.jsf? 
db=Tulcal). 

Transcripts from ML01 in the mycorrhizal roots were 
largely assigned to ribosomal proteins, enzymes that 
might be involved in plant cell wall degradation (1,4-beta- 
cellobiosidase A) and remodeling the fungal cell wall dur- 
ing growth and symbiosis (glycosyltransferase family 2 
protein, chitin deacetylase, UDP-glucose pyrophosphorylase, 
endo-beta-mannanase 4). We also identified numerous tran- 
scripts encoding different proteins associated with the nutri- 
ent transport and assimilation, such as H + -transporting 
ATPase, phosphate transporter (Unigene5906_AU), ammo- 
nium transporter (Unigene2660_All), amino acid and protein 
transporters, ZIP-like iron-zinc transporter, monosaccharide 
transporters (Unigene4644_All, Unigenel2394_AU) and 
MFS-type transporters (Unigenel6198_All, Unigene2162_All). 
ABC-type Fe 3+ -siderophore transporter and other transporters 
most probably implicated in the detoxification process 
were also detected in the mycorrhizal roots. The trans- 
lated amino acid sequence of Unigene5906_All shared 
59-60% identity to high-affinity Pi- transporters from 
G. versiforme (GvePT) and R. intraradices (RinPT) 
(Additional file 19: Figure Sll). The function of RinPT 
(also GinPT in references) has been suggested by Tisserant 
et al. [82] to be involved in stabilizing cooperation in the 
mycorrhizal symbiosis [82]. The transcript abundance of 
Unigene5906_All in CyEX23 was higher than that in CyEX21. 
So, this gene might be responsible for the better growth- 
promoting effect of MZ treatment than that of ML01 alone. 
Besides the transcripts in N and P assimilation pathway, iron- 
zinc transporter and ABC-type Fe 3+ -siderophore transporter 
co-existed in the mycorrhizal roots, which suggest an import- 
ant role of some metals in OM functionality. Our previous 
study showed that isolate CL01, Tulasnella-like Rhizocto- 
nia, formed typical mycorrhizae with orchids and en- 
hanced the N, P, Zn, Fe contents of mycorrhzial plants [8]. 
The ITS sequence of ML01 shared 84% similarity to that 
of isolate CL01 and both functioned similarly in improv- 
ing the growth and nutritional absorption of orchid plant- 
lets (our unpublished data). However, there is a need to 
compare the differences in gene expression between the 
intraradical and extraradical mycelium so as to elucidate 
the role of these transcripts in OM symbiosis. 



Conclusions 

Using a combination of RNA-seq and real-time RT-PCR, 
we characterized the root transcriptomic reprogramming 
of C. hybridum with respect to different beneficial fungi 
and identified specific changes in relative abundance for 
transcripts of genes in OM symbiosis. Genes involved in 
cell wall modification, reactive oxygen species detoxifica- 
tion, defense-related phytohormone and phosphate trans- 
port were co- induced in all the symbiotic interactions. 
The transcripts increased specifically by OM were related 
to signaling (LysM domain receptor-like kinase 3-like), 
cell wall degradation or reinforcement, protein metabol- 
ism and processing, defense (chitinase and mannose- 
specific lectin) and N, Fe transport. Among these DEGs, 
genes potentially functioned in C, N, P and Fe transport, 
hormone metabolism and signaling are of particular inter- 
est. Aside from these orchid transcripts, we also identified 
some putative fungal genes in symbiotic roots which were 
associated with plant cell wall degradation, remodeling the 
fungal cell wall and nutrient transport. This comprehen- 
sive dataset provides a basis for future research in explor- 
ing the molecular mechanisms underlying OM formation 
and functionality. 

Additional files 



Additional file 1: Table SI. Identities of the selected reference assemblies 
with Sanger-sequencing results of corresponding touch-down PCR products. 

Additional file 2: Table S2. Gene-specific primers used for qRT-PCR. 

Additional file 3: Figure SI. Typical co-cultivation condition of Cymbidium 
hybridum plantlets with different fungi at 15 dpi. Mock-inoculated control (A) 
and inoculated with isolate MI01 (B), ZH3A-3 (C) and ML01 + ZH3A-3 (D). 

Additional file 4: Figure S2. Root staining results of each treatment at 
different symbiotic time points. (A) Roots inoculated with E repens isolate 
ML01 at 6 days post-inoculation (dpi); (B) Roots inoculated with U. nana 
isolate ZH3A-3 at 6 dpi; (C) Roots inoculated with isolates ML01 and 
ZH3A-3 at 6 dpi; (D) Mock-inoculated roots at 6 dpi; (E) Roots inoculated 
with £ repens isolate ML01 at 10 dpi; (F) Roots inoculated with U. nana 
isolate ZH3A-3 at 10 dpi; (G) Roots inoculated with isolates ML01 and 
ZH3A-3 at 10 dpi; (H) Mock-inoculated roots at 10 dpi; (I) Roots inoculated 
with £ repens isolate ML01 at 15 dpi; (J) Roots inoculated with U. nana isolate 
ZH3A-3 at 15 dpi; (K) Roots inoculated with isolates ML01 and ZH3A-3 at 15 
dpi; (L) Mock-inoculated roots at 15 dpi; (M) Roots inoculated with £ repens 
isolate ML01 at 30 dpi; (N) Roots inoculated with U. nana isolate ZH3A-3 at 
30 dpi; (O) Roots inoculated with isolates ML01 and ZH3A-3 at 30 dpi; (P) 
Mock-inoculated roots at 30 dpi. 

Additional file 5. Cymbidium hybridum root transcriptome assembly- 
lllmina 90 bp reads. 

Additional file 6: Table S3. Highly abundant transcripts in CyTr-unigenes 
(Top 50). 

Additional file 7: Table S4. Pathway annotation of CyTr-unigenes. 

Additional file 8: Table S5. Statistics of DGE sequencing. 

Additional file 9: Table S6. Statistics of uniquely mapped reads for each 
sample. 

Additional file 10: Table S7-S9. Lists of differentially expressed genes 
in different C. hybridum- beneficial fungi interactions. 

Additional file 11: Table S10. Top 100 abundant genes in symbiotic 
roots with positive fold change in response to symbiotic fungi. Table S1 1 . 
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Top 100 abundant genes in non-symbiotic roots with negative fold change 
in response to symbiotic fungi. 

Additional file 12: Figure S3. Hierarchical clustering of DEGs from C. 
hybridum expressed in the roots in response to isolates ML01 (CyEX21), 
ZH3A-3 (CyEX22), ML01 +ZH3A-3 (CyEX23)or mock inoculation (CyEX20). 

Additional file 13: Lists of DEGS in 20 sub-clusters based on k-means. 

Additional file 14: Figure S4-S6. GO enrichment bar graphs of DEGs 
in the interactions between C. hybridum and ML01 (S4), ZH3A-3 (S5) or 
ML01 +ZH3A-3 (S6). 

Additional file 15: Table S12-S14. Lists of enriched pathways in different 
C. hybridum- beneficial fungi interactions. 

Additional file 16: Figure S7-S9. Scatterplot for top 20 enriched KEGG 
pathways of DEGs in the interactions between C hybridum and ML01 
(S7), ZH3A-3 (58) or ML01 +ZH3A-3 (59). 

Additional file 17: Figure S10. Comparison of expression profiles of 
random selected 26 genes by RNA-seq and qRT-PCR. 

Additional file 18: Table SI 5. Expression data for putative genes from 
£ repens isolate ML01 in symbiotic roots. Table S16. Expression data for 
putative genes from U. nana isolate ZH3A-3 in symbiotic roots. 

Additional file 19: Figure S11. Phylogenetic tree for the amino acid 
sequences of symbiosis- associated phosphate transporters from plant 
and fungi.The dendrogram was generated by Mega 4.0 software using 
ClustlW for the alignment and the neighbor-joining method for the 
construction of phytogeny. The plant Pi-transporters accession numbers 
could be obtained in [49]. Fungal phosphate transporters accession 
numbers: Glomus versiforme: GvePT, AAC49132.1; GvePT2, gi|74654973; 
Rhizophagus intraradices: RinPT, AAL37552.1 ; Yeast: YEPH084, gi| 
1346710: Sesbania rostrata: SerPT, CAC28219.1; Piriformospora indica: 
PiPTl, ABI93950.1; PiPT2, CCA74169.1; PiPT3, CCA74178.1; PiPT4, 
CCA76379.1 . The nucleotide sequences of putative phosphate 
transporters in our study were available in Additional file 6: Table S3. 
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